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Pattern formation and evolution in unsynchronizable complex networks are investigated. Due to the asym- 
metric topology, the synchronous patterns formed in complex networks are irregular and nonstationary. For 
coupling strength immediately out of the synchronizable region, the typical phenomenon is the on-off inter- 
mittency of the system dynamics. The patterns appeared in this process are signatured by the coexistence of a 
giant cluster, which comprises most of the nodes, and a few number of small clusters. The pattern evolution 
is characterized by the giant cluster irregularly absorbs or emits the small clusters. As the coupling strength 
leaves away from the synchronization bifurcation point, the giant cluster is gradually dissolved into a number of 
small clusters, and the system dynamics is characterized by the integration and separation of the small clusters. 
Dynamical mechanisms and statistical properties of the nonstationary pattern evolution are analyzed and con- 
ducted, and some scalings are newly revealed. Remarkably, it is found that the few active nodes, which escape 
from the giant cluster with a high frequency, are independent of the coupling strength while are sensitive to the 
bifurcation types. We hope our findings about nonstationary pattern could give additional understandings to 
the dynamics of complex systems and have implications to some real problems where systems maintain their 
normal functions only in the unsynchronizable state. 



PACS numbers: 89.75.-k, 05.45.Xt 

I. INTRODUCTION 

Synchronization of complex networks has aroused many in- 
terest in nonlinear science since the discoveries of the small- 
world and scale-free properties in many real and man-made 
systems (U, 0]. In this study, one important issue is to ex- 
plore the inter-dependent relationship between the collective 
behaviors of the complex systems and their underlying topolo- 
gies. In particular, many efforts have been paid to the con- 
struction of optimal networks, and a number of factors which 
have important affections to the synchronizability of complex 
networks have been gradually disclosed. Now it is known that 
random networks, due to their small average distances, are 
generally more synchronizable than regular networks |H,3]; 
and scale-free networks, with weighted and asymmetric cou- 
plings, can be more synchronizable than homogeneous net- 
works H,^. In these studies, the standard method employed 
for synchronization analysis is the master stability function 
(MSF), where the network synchronizability is estimated by 
an eigenratio calculated from the coupling matrix, and system 
which has a smaller eigenratio is believed to be more synchro- 
nizable than that of larger eigenratio [7]. Inspired by this, to 
improve the network synchronizability, the only task seems to 
be upgrading the coupling matrix so as to decrease the eigen- 
ratio, either by changing the network topology (4j] or by ad- 
justing the coupling scheme |5|,|g]- 

The MSF method, while bringing great convenience to the 
analysis, overlooks the temporal, local property of the sys- 
tem and reflects only partial information about the system 
dynamics. Specifically, from MSF we only know ultimately 
whether the network is globally synchronizable or unsynchro- 
nizable, but do not know how the global synchronization is 



reached if the network is synchronizable, or what's the pat- 
tern and how it evolves if the network is unsynchronizable. 
These evolution details, or the transient behavior in system 
development, contain rich information about the system dy- 
namics and may give additional insights to the organization 
of complex systems. For instance, the recent studies about 
synchronization transition have shown that, in the unsynchro- 
nizable states, heterogeneous networks are more synchroniz- 
able (have a higher degree of coherence) than homogeneous 
networks at small couplings, while at larger couplings the op- 
posite happens [8]. This crossover phenomenon of network 
synchronizability are difficult to understood if we only look 
at the final state of the system, but are straightforward if we 
look at the transient behaviors of their evolutions |8]. Besides 
revealing the synchronization mechanisms, the transient be- 
havior of network synchronization can also be used to detect 
the topological scales and hierarchical structures in the real 
systems, e.g., the detection of cluster structures in social and 
biological networks IHQii. However, despite of its theoreti- 
cal and practical significance, the study of transient dynamics 
of complex networks is still at its infancy and many questions 
remain open, say for example, the pattern evolution of unsyn- 
chronizable complex networks. 

Pattern formation in unsynchronizable but near- 
synchronization networks has been an important issue 
in studying the collective behavior of regular networks 
ifTll Wm . By setting the coupling strength nearby the 
synchronization bifurcation point, the system state shares 
both the dynamical properties of the synchronizable and 
unsynchronizable states: a state of high coherence but is not 
synchronized. The bifacial dynamical property makes this 
state a natural choice in investigating the transition process of 
networks synchronization. Previously studies about regular 



2 



networks, say for example the lattices flUl . have shown that, 
when the coupling strength is slightly out of the synchroniz- 
able region, although global synchronization is unreachable, 
nodes are still synchronized in a partial sense. That is, nodes 
are self-organized into a number of synchronous clusters. 
The distribution of these clusters, also called the synchronous 
pattern, is determined by a set of factors such as the coupling 
strength, the system size and the coupling scheme. As the 
coupling strength leaves away from the bifurcation point, the 
pattern structure becomes more and more complicated and the 
system coherence will be decreased, and finally reaches the 
turbulence state. It is worthy of note that the patterns arisen 
in regular networks have two common properties: spatially 
symmetric and temporally stationary. More specifically, the 
contents of each cluster are fixed and the clusters are of 
translation symmetry in space. For this reason, we say that 
the patterns formed in regular networks are symmetric and 
stationary. These two properties, as have been discussed in 
the previous studies 111 ll Il2ll . are rooted in the symmetric 
topology of the regular networks. This makes it interesting 
to ask the following question: how about the patterns in 
unsynchronizable complex networks? 

Different to the regular networks, in complex networks we 
are not able to find any symmetry from their topologies. The 
asymmetric topology, according to the pattern analysis devel- 
oped in studying regular networks ll ill , will induce two sig- 
nificant changes to the patterns: 1) the synchronous clusters, if 
they exist, will be asymmetric; and 2) all the possible patterns, 
including the one of global synchronization, are linearly un- 
stable under small perturbations. In other words, the patterns 
formed in complex networks are expected to be asymmetric 
and nonstationary. Our mission of this paper is just to under- 
stand and characterize the nonstationary patterns arisen in the 
development of complex networks. Specifically, we are trying 
to investigate the following questions: 1) is there any pattern 
arises during the system evolution? 2) the pattern is stationary 
or nonstationary? if nonstationary, how is it evolving and how 
is it reflected from the system dynamics? 3) What happens 
to the pattern properties during the transition of network syn- 
chronization? and 4) How the coupling strength and bifurca- 
tion type affect the pattern properties? By investigating these 
dynamical and statistical properties, we wish to have a global 
understanding to the dynamics of unsynchronizable complex 
networks. 

Our main findings are: 1) for coupling strength immedi- 
ately outside of the synchronizable region, the system dynam- 
ics undergoes the process of on-off intermittency. That is, 
most of the time the system stays on the global synchroniza- 
tion state (the "off state) but, once in a while, it develops 
into a breaking state (the "on" state) which is composed by 
a giant cluster and a few number of small clusters (hereafter 
we call it the giant-cluster state). As the system develops, 
the giant cluster changes its shape by absorbing or emitting 
the small clusters, leading to the "off or "on" states, respec- 
tively; 2) the few active nodes which escape from the giant 
cluster with the high frequencies are coupling-strength inde- 
pendent but are bifurcation-type dependent. That is, in the 
neighboring region of a fixed bifurcation point, the locations 



of these active nodes do not change with the coupling strength; 
if we change the coupling strength from nearby another bi- 
furcation point (the two bifurcation points will be explained 
later), their locations will be totally changed; 3) as coupling 
strength leaves away from the bifurcation point, the giant clus- 
ter is gradually dissolved and more small clusters are gener- 
ated from it. Eventually, the giant cluster disappears and the 
pattern is composed by only the small clusters (hereafter we 
call it the scattering-cluster state). During the course of sys- 
tem evolution, each small cluster may either increase its size 
by integrating with other small clusters or decrease its size 
by breaking to even small clusters, but it can never reach to 
the global synchronization state; 4) besides the giant cluster, 
the giant- and scattering-cluster states are also distinct in their 
small clusters. For giant-cluster state the size of the small 
clusters follows a power-law distribution, while for scattering- 
cluster state it follows a Gaussian distribution. 

The rest of the paper is going to be arranged as follows. In 
Sec. II we will give our model of coupled map network and, 
based on the method of MSF, point out the two bifurcation 
points and the transition areas that we are going to study with. 
In Sec. Ill we will employ the method of finite-time Lyapunov 
exponent to predict and describe the intermittent system dy- 
namics in the bifurcation regions. Direct simulations about 
on-off intermittency will be presented in Sec. IV. By intro- 
ducing the method of temporal phase synchronization, in Sec. 
V we will investigate in detail the dynamical and statistical 
properties of the nonstationary pattern. Meanwhile, proper- 
ties of the giant- and scattering states will be compared and 
the transition between the two states will be conducted. In 
Sec. VI we will discuss the phenomenon of active nodes and 
investigate their dependence to the network properties. Dis- 
cussions and conclusions about pattern evolution in complex 
networks will be presented in Sec. VII. 

II. COUPLED MAP NETWORKS AND THE TWO 
BIFURCATION POINTS 

Our model of coupled map network is of the following form 

x,(f + 1) = F(xi(t)) - e G id U [f (x, (t))} . (1) 

3 

where x,(t + 1) = F(xj(t)) is a d-dimensional map repre- 
senting the local dynamics on node i, e is a global coupling 
parameter, G is Laplacian matrix representing the couplings, 
and H is a coupling function. To facilitate our analysis, we 
adopt the following coupling scheme Hit]: 

Aij k j 

G ^ = A (2) 

for j ^ i and d.i = 1, with fcj the degree of node i and 
A the adjacent matrix of the network: Ajj = 1 if node i and 
j are connected and Aij = otherwise. In comparison with 
the traditional coupling schemes, one important advantage we 
benefit from this coupling scheme is that the synchronizabil- 
ity of the network, i.e. the eigenratio of the coupling matrix 
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described in Eq. Q, can be easily adjusted by the parameter 
j3, while the network topology is kept unchanged. This advan- 
tage brings many convenience in network selection since for 
a given network topology, even though it is unsynchronizable 
under the traditional schemes, can now be synchronizable by 
adjusting j3 in Eq. Q. This convenience is of particular im- 
portance when our studies of network dynamics are focused 
on the bifurcation regions, where the network synchronizabil- 
ity should be deliberately arranged in order to demonstrate 
both the two types of bifurcations. We note that the adop- 
tion of Eq. J2) is only for the purpose of convenient analy- 
sis, the findings we are going to report are general and can 
also be observed by other coupling schemes given the net- 
work is properly prepared. In practice, we use logistic map 
F(x) = 4x(l — x) as the local dynamics and adopt H(x) = x 
as the coupling function. 

We first locate the two bifurcation points of global syn- 
chronization. The linear stability of the global synchroniza- 
tion state {xi(t) — s(t),Vi} is determined by the correspond- 
ing variational equations, which can be diagonalized into N 
blocks of form 

y(t + l) = [T>F(s)+aT>n(s)]y(t), (3) 

with DF(s) and DH(s) the Jacobian matrices of the corre- 
sponding vector functions evaluated at s(t), and y represents 
the different modes that are transverse to the synchronous 
manifold s(t). We have a(i) = e\ for the ith block, i — 
1, 2, N, and Ai = < A2 < ... < \n are the eigenvalues 
of matrix G. The largest Lyapunov exponent A(er) of Eq. Q, 
known as the master stability function (MSF) [7], determines 
the linear stability of the synchronous manifold s(t). In par- 
ticular, the synchronous manifold is stable if only A(eAj) < 
for each i = 2, N, The set of Lyapunov exponents A(eXi) 
govern the stability of the synchronous manifold in the trans- 
verse spaces, and a positive value of A(eAj) represents the 
loss of the stability in the transverse space of mode i. It was 
found that for a large class of chaotic systems, A(er) < 
is only fulfilled within a limit range in the parameter space 
a G (<7i, 02). This indicates that, to make the global synchro- 
nization state linearly stable, all the eigenvalues Ai should be 
contained within range (01,02), i.e., < 02/01. For 

the logistic map employed here, it is not difficult to prove that 
ai = 0.5 and 02 = 1.5. Therefore, to achieve global syn- 
chronization, the coupling matrix G should be designed with 
eigenratio R = Ajy/A 2 < 02/01 = 3 = R c . 

Besides the condition of R < R c , to guarantee the synchro- 
nization, we also need to set the coupling strength in a proper 
way: either small or large couplings may deteriorate the syn- 
chronization. If e < s\ = 01/A2, the couplings are too weak 
to restrict the node trajectories to the synchronous manifold; 
while if e > £2 = 02/ Ajv. the couplings will be too strong 
and actually act as large perturbations to the synchronization 
manifold. Therefore, to achieve the global synchronization, 
we also require e\ < e < £2- The two critical couplings 
E\ and £2, which are named as the long-wave (LW) 01411 and 
short-wave (SW) bifurcations [15] respectively in the studies 
of regular networks, thus stand as the boundaries of the syn- 
chronizable region. Our studies about network synchroniza- 
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FIG. 1: For scale-free network of N — 1000 nodes and of average 
degree < k >= 8, a schematic plot on the generation of the two 
bifurcation points as a function of the coupling strength. The long- 
wave bifurcation occurs at about ei « 0.83488 which is determined 
by the condition e\2 = 01 = 0.5 (the lower line). The short-wave 
bifurcation occurs at about £2 ~ 0.95127 which is determined by 
the condition e\n = 02 — 1.5 (the upper line). 

tion will be focused on the neighboring regions of the two 
bifurcation points, i.e., the region of e < E\ or e > £2. 

By the standard BA growth model [ 1], we construct a scale- 
free network of 10 3 nodes and of average degree (k) = 8. By 
setting [3 — 2.5 in Eq. Q, we have A2 ~ 0.6 and A^v ~ 1-58. 
Because of R = ~ 2.6 < R c , the network is glob- 

ally synchronizable. Also, because of A2 > <J\ and Xn > 0-2, 
both the two bifurcations can be realized by adjusting the cou- 
pling strength within range e G (0, 1). In specific, when 
£ < £1 ~ 0.835, we have £A2 < cr-y and £A^ < 02, the 
synchronous manifold loses its stability at the lower boundary 
of the synchronizable region and LW bifurcation occurs; and 
when £ > £2 ~ 0.95, we have £A2 > o\ and eXn > 02, the 
synchronous manifold loses its stability at the upper bound- 
ary of the synchronizable region and SW bifurcation occurs 
[Fig. 1]. In the following we will fix the network topology 
and the parameter (3, while generating the various patterns by 
changing the coupling strength e nearby the two bifurcation 
points. 



III. FINITE-TIME LYAPUNOV EXPONENT 

Before direct simulations, we first give a qualitative de- 
scription (prediction) on the possible system dynamics in bi- 
furcation regions. To concrete our analysis, in the following 
we will only discuss the situation of SW bifurcation (e < £1), 
while noting that the same phenomena can be found at the 
LW bifurcation as well (e > £2) . In preparing the unsynchro- 
nizable states, we only let A(A2) be slightly puncturing into 
the unstable region, while keeping all the other exponents still 
staying in the stable region, i.e., A(A2) > and A(A,-) < 
for i — 3, ...N. With this setting, the synchronous manifold 
is only desynchronized in the transverse space of mode 2. As 
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such, the system possesses only two positive Lyapunov expo- 
nents, one is A(Aq) which is associated to the synchronous 
manifold itself and another one is A(A2). Noticing that A(A) 
are asymptotic averages, and, as so, they account only for the 
global stability properties, but do not warrant the possible co- 
herent sets arising in the system evolutions. These coherent 
sets, for regular networks, refer to the stationary, symmetric 
patterns to which the system finally develops. While for com- 
plex networks, these sets can be the temporal, irregular clus- 
ters emerged in the process of system evolution. 

In the region of e < s\, although global synchronization is 
unreachable, the system may still keep with the high coher- 
ence due to the existence of the synchronous clusters. Espe- 
cially, there could be some moments at which all the trajecto- 
ries are restrained to a small region in the phase space, very 
close to the situation of global synchronization. This vary- 
ing system coherence, however, can not be reflected from the 
asymptotic value A(A). To characterize the variation, we need 
to employ some new quantities which are able to capture the 
temporal behavior of system. One of such quantities is the 
finite-time Lyapunov exponent (FLE), a technique developed 
in studying chaos transition in nonlinear science 11611 . In stead 
of asymptotic average, FLE measures the diverging rate of 
nearby trajectories only in a finite time interval T. 

1 lT 

A i = - l*DH(«(t)). (4) 

t=(i-l)T 

As our studies are focused on the situation of one-mode 
desynchronization, the stability of the synchronous manifold 
and the temporal behavior that it displays are therefore ex- 
pected to be more reflected from the variation of Aa », the 
FLE that associates with mode 2. With e — 0.83 and 
T = 100, we plot in Fig. 2 the time evolution of A2,,. It 
is found that, although with a positive asymptotic value about 
(A2.1) ~ 6 x 10~ 3 , the instant value of Ag.t penetrates to the 
negative region at a high frequency. According to the different 
signs of Aa,i, the system evolution is divided into two types of 
intervals: the divergent interval and the contractive interval. 
In the divergent intervals we have A2,, > and the system 
dynamics is temporarily dominated by the divergence of the 
node trajectories from the synchronous manifold; while in the 
contractive intervals we have A2,j < and the system dynam- 
ics is temporarily dominated by the convergence of the node 
trajectories to the synchronous manifold. 

The variation of A2 j, reflected on the process of pattern 
evolution, characterizes the travelling property of the sys- 
tem dynamics among the neighboring regions of two different 
kinds of states: the desynchronization state and the synchro- 
nization state. In Fig. 2, the minimum value of Ag.j is about 
—0.07, during this contractive interval the node trajectories 
will converge to the synchronous manifold by an amount of 
e mmA 2 ,iT w g -7 w io~ 3 on average. Assuming that before 
entering this interval the average distance between the node 
trajectories is A (for logistic map we always have A < 1), 
then at the end of this interval the average distance is de- 
creased to A x 10 -3 , a small value which is usually over- 
shadowed by noise in practice. Due to this small distance, 
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FIG. 2: For e = 0.83 in Fig. 1, the time evolution of the finite- 
time Lyapunov exponent Aa,« calculated on intervals of length T = 
100. It is observed that, while having the positive asymptotic value 
(Aa,i) > 0, the temporal value of Aa.j is penetrating into the negative 
region frequently. 

the system can be practically regarded as already reached the 
synchronization state. On the other hand, if the system enters 
a divergent interval, the node trajectories will diverge from 
each other and, at the end of this interval, their average dis- 
tance will be increased by an order of 10 3 . This large distance 
will deteriorate the ordered trajectories (or the high coherence 
of the system dynamics) that achieved during the contractive 
intervals, and leading to the incoherent, breaking state. The 
pattern of the breaking state, however, is not unique. Depend- 
ing on the initial conditions and the divergence intervals, the 
pattern may assume the different configurations. Therefore, 
based on the observations of Aa,i [Fig. 2] the dynamics of 
unsynchronizable networks can be intuitively understood as 
an intermittent travelling among the synchronization state and 
the different kinds of desynchronization states. 



IV. ON-OFF INTERMITTENCY DESCRIBED BY 
COMPLETE SYNCHRONIZATION 

We now investigate the system dynamics by direct simula- 
tions. To implement, we first prepare the system to be staying 
on the synchronization state. This can be achieved by adopt- 
ing a large coupling strength from the synchronizable region, 
i.e. £1 < e < £2- After synchronization is achieved, we then 
decrease £ to a value slightly below the bifurcation point E\ 
and, in the meantime, an instant small perturbation is added 
on each node. In practice, we take i.i.d (independent iden- 
tically distributed) noise of strength 1 x 10~ 5 as the pertur- 
bations. After this, we release the system and let it develop 
according to Eq. (Q]). Since e < £1, the synchronization state 
is unstable and, triggered by the noises, the node trajectories 
begin to diverge from each other. The divergent trajectories, 
however, will frequently visit the neighborhood of the syn- 
chronous manifold, especially during those contractive inter- 
vals of small A2.i [Fig. 2]. The intermittent system dynamics 
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is plotted in Fig. 3(a), where the average trajectory distance 
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N 



x is plotted as a function of time. As 



AX 

we have predicted from LLE, the system dynamics indeed un- 
dergoes an intermittent process. To characterize the intermit- 
tency, we plot in Fig. 3(b) the laminar-phase distribution of 
the AX sequence plotted in Fig. 3(a). It is found that the lam- 
inar length t (the time interval between two adjacent bursts 
of amplitude AX(t) > 10 -3 ) and the probability p(r) for it 
to appear follow a power-law scaling p(r) ~ t~ 7 . The fitted 
exponent is about 7 w — 1.5 ± 0.05, with a fat tail at large r 
due to the finite simulating time. 

In chaos theory, intermittent process of laminar-phase expo- 
nent —3/2 is classified as the "on-off ' intermittency, a typical 
phenomenon observed in dynamical systems with a symmet- 
ric invariant set 117!]. On-off intermittency is also reported 
in chaos synchronization of regular networks, where the in- 
variant set refers to the synchronous manifold, and the "off 
state refers to the long stretches that the system dynamics is 
staying nearby the synchronous manifold and the "on" state 
refers to the short bursts that the system dynamics is staying 
away from the synchronous manifold. Therefore, in terms of 
laminar-phase distribution, the intermittency we have found 
in complex networks [Fig. 3] has no difference to the that of 
the regular networks, despite of the drastic difference between 
their topologies. We have also investigated the transition be- 
havior of the averaged distance (AX(t)) nearby the bifurca- 
tion points. As shown in Fig. 3(c), a linear relation between 
(AX(t)) and e is found in the region of e < E\. This linear 
transition of the system performance, again, is consistent with 
the transition of regular networks 111811 . Therefore, in terms of 
complete synchronization, the on-off intermittency we have 
found in complex networks has no difference to that of the 
regular networks. 



V. PATTERN EVOLUTION IN COMPLEX NETWORKS 

To reveal the unique properties of the system dynamics that 
induced by the complex topology, we go on to investigate 
the pattern formation of unsynchronizable networks by the 
method of temporal phase synchronization (TPS). 



A. Temporal phase synchronization 

TPS is defined as follows. Let Xi(t) be the time sequence 
recorded on node i, we first transform it into a symbolic se- 
quence 8i(t) according to the following equations 



if Xi(t) < 0.5, 
if Xi{t) > 0.5. 



(5) 



Then we divide 6i(t) into short segments of the equal length 
n. Regarding each segment as an new element, we there- 
fore have transformed the long, variable sequence Xi (t) into 
a short, symbolic sequence 0;(i'). If at moment t' we have 
Qi(t') = <dj(t'), then we say that TPS is achieved between 
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FIG. 3: (Color online) The on-off intermittency of the system dy- 
namics nearby the LW bifurcation at e = 0.83. (a) The time evolu- 
tion of the average trajectory distance AX. (b) The laminar-phase 
distribution of AX, which follows a power-law scaling with the fit- 
ted exponent around 3/2. (c) The transition behavior of the average 
distance (AX) nearby the LW bifurcation point e%, where a linear 
relation is found between the two quantities. 



the nodes i and j. The collection of nodes which have the 
same value of Q at moment t' are defined as a temporarily 
synchronous cluster, and all the synchronous clusters consti- 
tute the temporarily pattern of the system. During the course 
of system evolution, the clusters will change their shapes and 
contents and the pattern will change its configuration. 

In comparison with the method of complete synchroniza- 
tion, the advantage we benefit from TPS is obvious: it makes 
the synchronous pattern detectable. With complete synchro- 
nization, it is almost impossible for two nodes to have ex- 
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actly the same variable at the same time. Despite the fact that 
at some moments the system has already reached the high- 
coherence states (formed during those contractive intervals in 
Fig. 3(a)), with complete synchronization we are not able to 
distinguish these states from those low-coherence ones quanti- 
tatively (formed during those divergent intervals in Fig. 3(a)). 
(A remedy to this difficulty seems to define the clusters by 
the method of threshold truncation, i.e., nodes are regarded 
as synchronized if the distance between their trajectories is 
smaller than some small value. However, this definition of 
synchronization will induce the problem of cluster idenfica- 
tion, as the same state may generate different patterns if we 
choose the different reference nodes.) On the contrary, TPS 
focuses on the loose match (phase synchronization) between 
the node variables over a period of time. By requiring an ex- 
act match of the discrete variable 6, the synchronous pattern 
is uniquely defined; while by requiring the match of the long 
sequences of 8, the "synchronous" nodes are guaranteed with 
a strong coherence. 



B. Pattern evolution of the giant-cluster state 

With the same set of parameters as in Fig. 3(a), by the 
method of TPS we plot in Fig. 4 the time evolutions of two ba- 
sic quantities of pattern evolution: the number of synchronous 
clusters n c and the size of the largest cluster L max . It is 
found that, similar to the phenomenon in complete synchro- 
nization [Fig. 3(a)], on-off intermittency is also found in the 
TPS quantities n c and L max . In Fig. 4(a) it is shown that most 
of the time the system is broken into only a few number of 
clusters, n c = 2 or 3, while occasionally it is broken into a 
quite large number of clusters, 10 < n c < 50, or united to the 
synchronization state, n c = 1. The intermittent pattern evo- 
lution is also reflected on the sequence of £ max [Fig. 4(b)], 
where most of the time the size of the giant cluster is about 
imax ~ N, while occasionally it decreases to some small 
values of L max < N/2. As we have discussed previously, 
the main advantage we benefit from TPS is in identifying the 
clusters. This advantage is clearly shown in Figs. 4(a) and 
(b), where for any time instant the two quantities n c and Li 
are uniquely defined. Besides cluster identification, we also 
benefit from TPS in quantifying the synchronization degrees. 
In specific, the different coherence states shown in Fig. 3(a) 
now can be clearly quantified: high coherence states are those 
of smaller n c or larger £ max . Specially, the synchronization 
state now is unambiguously defined as the moments of n c = 1 
in Fig. 4(a) or, equally, the moments of L max = N in Fig. 
4(b). 

We go on to investigate the pattern evolution by statistical 
analysis. The first statistic we are interested is the laminar- 
phase distribution of the synchronization state, i.e. the time 
intervals that n c = 1 in Fig. 4(a) or L max = N in Fig. 4(b). 
In its original definition, laminar phase refers to the time inter- 
vals t that all node trajectories stays within a small distance 
from the synchronous manifold, therefore the actual value of 
r is varying with the predefined threshold distance. This un- 
certainty is overcome in TPS. As shown in Fig. 4(a), in TPS 
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FIG. 4: For the same set of parameters as in Fig. 3(a). The time 
evolutions of the TPS quantities, (a) The number of the synchronous 
clusters n c and (b) the size of the giant cluster Lmax- The synchro- 
nization state is defined as the moments n c = 1 in (a) or L max = N 
in (b). 



the "off state refers to the moments of n c = 1 specifically. 
The laminar-phase distribution of n c is plotted in Fig. 5(a). In 
consistency with the distribution of complete synchronization 
[Fig. 3], the laminar-phase distribution of n c also follows a 
power-law scaling and has the same exponent 7 w — 1.5±0.1. 
Therefore the use of TPS, while bringing convenience to the 
pattern analysis, still capture the basic properties of the on-off 
intermittency. The second statistic we are interested is the size 
distribution of the largest cluster, an important indicator for 
the coherence degree of the system. For the L max sequence 
plotted in Fig. 4(b), in Fig. 5(b) we plot the size distribution 
of L max . It is seen that the probability of finding large clus- 
ter I/ max « N is much higher than that of small cluster of 
L max < 500. In particular, the probability for finding clusters 
of L max > 990 is about 20 percent and for L max > 990 it 
is about 70 percent. Therefore, in the region of s < £1, the 
distinct feature of the system patterns is the existence of a gi- 
ant cluster. Due to this special feature, we call these states the 
giant-cluster state. 

Besides the giant cluster, we are also interested in the prop- 
erties of the small clusters. We plot in Fig. 5(c) the distribu- 
tion of n c and in Fig. 5(d) the size distribution of the small 
clusters L, that surround the giant cluster in the pattern. As 
shown in Fig. 5(c), the distribution of n c follows a power-law 
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FIG. 5: (Color online) Statistical properties of the on-off intermit- 
tency plotted in Fig. 4. (a) The power-law scaling of the laminar- 
phase distribution of n c . The fitted slope is about —2.3 ± 0.05. (b) 
The size distribution of the size of the giant cluster, (c) The power- 
law distribution of the number of small clusters n c . The fitted slope 
is about —3 ± 0.1. (d) The power-law scaling on the size distribution 
of the small clusters. The fitted slope is about —1.2 ± 0.01. 



scaling with the fixed exponent is about 7 « —3 ± 0.05. The 
heterogeneous distribution of n c indicates that, in the giant- 
cluster state, the system is usually broken into only a few num- 
ber of clusters. An interesting finding exists in the size distri- 
bution of the small clusters. As shown in Fig. 5(d), in range 
L{ £ [1, N/2] a power-law scaling is found between Pn and 
Li, with the fitted exponent is about 7 sw —1.1 ± 0.05. The 
distribution of L, confirms the finding of Fig. 5(c) that the 
small clusters which join or separate from the giant cluster are 
usually of small size. 

Combining the findings of Fig. 4 and Fig. 5, the picture 
of pattern evolution in the bifurcation region e < E\ now be- 
comes clear. Generally speaking, the evolution can be divided 
into two opposite dynamical processes happening around the 
giant cluster: the separation and integration of the small clus- 
ters. During the separation process, the small clusters are es- 
caped from the giant cluster, which weakens the dominant role 
of the giant cluster and makes the pattern complicated. How- 
ever, the separated clusters occupy only a small proportion of 
the nodes [Fig. 5(c)], the majority nodes are still attached to 
the giant cluster, which sustains the synchronization skeleton 
and keeps the system on the high coherence states. At some 
rare moments the giant cluster may disappears, and the pat- 
tern is composed by only small clusters of Li < N/2. At 
these moments, the synchronization skeleton is broken, the 
pattern becomes even complicated and the system coherence 
reaches its minimum. In contrast, during the process of cluster 
integration, the giant cluster will increase it size by attracting 
the small clusters , and gradually towards the state of global 
synchronization. It should be noticed that the separation and 
integration processes are uneven and are typically occurring 
at the same time. For instance, during the separation process, 
while the system evolution is dominated by the separation of 



new small clusters from the giant cluster, there could be some 
small clusters rejoin to the giant cluster. 



C. Pattern evolution of the scattering-cluster state 

As we further decrease the coupling strength from e%, the 
picture of pattern evolution will be totally changed. With 
e = 0.79, we plot in Fig. 6 the same statistics as in Fig. 5. 
The first observation is the loss of the global synchronization 
state, as can be found from the time variation of n c plotted in 
Fig. 6(a). The loss of global synchronization becomes even 
clear if we compare Fig. 6(a) with Fig. 4(a): in Fig. 6(a), 
except the moment at t — 0, the system can never reach the 
synchronization state at n c = 1 and very often it is broken into 
a large number of small clusters at about n c ~ 10 2 . The fact 
that the pattern is decomposed into a large number of small 
clusters is also manifested by the distribution of n c , as plotted 
in Fig. 6(b). Instead of the power-law distribution found in the 
giant-cluster state, in the scattering-cluster state n c follows a 
Gaussian distribution [Fig. 6(b)]. As e further decreases from 
El, the mean value of n c will shift to the larger values, as 
indicated by the e = 0.78 curve plotted in Fig. 6(b). The sec- 
ond observation is the disappearance of the giant cluster. As 
shown in Fig. 6(c), the size distribution of the largest cluster 
also follows a Gaussian distribution, with its mean value lo- 
cates at (L mra ) < N/2. The distribution of Fig. 6(c) is very 
different to that of Fig. 5(b), where in Fig. 5(b) the largest (gi- 
ant) cluster has size i max N in most of the time. As e de- 
creases, the mean value of the largest cluster (L max ) will shift 
to small values and the variance of L max will be decreased, as 
indicated by the e = 0.78 curve plotted in Fig. 6(c). Similar to 
plot of Fig. 5(d), we have also investigated the distribution of 
Li, the sizes for all the small clusters appeared in the system 
evolution [Fig. 6(d)]. It is found that the distribution of Li fol- 
lows a power-law distribution for Li < N/2, while having an 
exponential tail for Li > N/2. Numerically we find that the 
exponent of the power-law section, i.e. in range Li £ [1,200), 
is about —2 ± 0.05, while the fitted exponent for the exponen- 
tial section is about —4.5 x 10 -3 ±2 x 10~ 5 . These two 
exponents, however, are changing with e. As e decreases, the 
two exponents will shift to some small values. 

Combining Fig. 5 and Fig. 6, we are able to outline the 
transition process of network synchronization nearby the bi- 
furcation points, i.e., the transition from the giant-cluster state 
to the scattering-cluster state as e leaves away from e\. In 
the region of e < E\, the pattern is composed by a giant 
cluster and a few number of small clusters, i.e. the giant- 
cluster state. As e decreases from e\ gradually, more and 
more small clusters will be emitted out from the giant clus- 
ter and, as a consequence, both the size of the giant cluster 
and the fraction of synchronization time will be decreased. 
Then, at about e c « 0.832, the giant cluster disappears and 
the pattern of the system is composed by several larger clus- 
ters, of size L max < N/2, together with many small clusters 
of heterogenous size distribution, i.e. the scattering-cluster 
state. After that, as e decreases from e c , the clusters shrink 
their size by breaking into even small clusters, and the pat- 
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FIG. 6: (Color online) The dynamical and statistical properties of 
pattern evolution for e = 0.79. (a) The time evolution of n c . (b) The 
Gaussian distribution of the number of the small clusters n c . (b) The 
Gaussian distribution of the size of the largest cluster I/ m ax- (d) The 
two-segment distribution on the size of the small clusters Li. In the 
region of Li < 200, Li follows a power-law distribution with fitted 
exponent is about —2 ± 0.05; while for Li > 200, the distribution 
is exponential with the fitted exponent is about —4.5 x 10~ 3 ± 2 x 
10 . As e further decreases from ei, the largest cluster becomes 
even smaller and more small clusters are emitted out from it. As 
illustrated by the e = 0.78 curves plotted in (b), (c) and (d). 



tern becomes even complicated. The detail transition from the 
giant-cluster state to the scattering-cluster state is presented 
in Fig. 7, where the average number of clusters that the sys- 
tem is broken into (n c ), Fig. 7(a), and the average size of 
the largest cluster (L max ), Fig. 7(b), are plotted as a func- 
tion of the coupling strength in the LW bifurcation region. 
The transition is found to be smooth and steady, just as we 
have expected. Besides the giant cluster, another difference 
between the giant-cluster and scattering-cluster states exists 
in their pattern evolutions. In the giant-cluster state, while 
the configuration of the giant cluster is continuously updated 
by emitting or absorbing the small clusters, its main contents 
are stable and do not change with time. In contrast, in the 
scattering-cluster state the small clusters integrate with or sep- 
arate from each other in a random fashion. Although occa- 
sionally there could be some large clusters show up in the pat- 
tern of the scattering-cluster state [Fig. 6(d)], these "large" 
clusters, however, are very fragile and will break into small 
clusters again in a short time. This quick-dissolving prop- 
erty stops the scattering-cluster state from having a high co- 
herence. 



VI. CHARACTERIZING THE ACTIVE NODES 

In the giant-cluster state, most of the nodes are organized 
into the giant cluster while few nodes, either in forms of small 
group or isolated node, are separating from or joining to the 
giant cluster with a high frequency. These active nodes, al- 
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FIG. 7: The transition process of the network synchronization nearby 
the LW bifurcation point ei. (a) The average number of clusters that 
the system is broken into as a function of coupling strength, (b) The 
average size of the largest cluster as a function of coupling strength. 
Each date is an averaged result over 10 8 time steps. 



though are few in amount, play an important role in network 
synchronization. Clearly, a proper characterization of these 
nodes will deepen our understandings on the system dynam- 
ics and give indications to the improvement of network perfor- 
mance. For instance, to improve the synchronizability of the 
system, we may either remove the few most active nodes from 
the network, or update their coupling strengths specifically. 

In characterizing the active nodes, the following properties 
are of general interest: 1) what's the dependence of the node 
activity on the network topology? can we characterize these 
nodes by the known network properties such as node degree 
orbetweenness? 2) are their locations sensitive to the coupling 
strength? and 3) what's the effect of bifurcation type on their 
locations? In the following we will explore these questions by 
numerical simulations. 

We first try to characterize the active nodes by their topo- 
logical properties. For the giant-cluster state described in Fig. 
4, we plot in Fig. 8(a) the probability p u \ that each node stays 
in the giant cluster. While the majority nodes stay in the gi- 
ant cluster with a high probability p u i w 1, few nodes are 
of unusually small probabilities: 1 percent of the nodes have 
Pui < 0.8. One important observation of Fig. 8(a) is that 
the locations of the active nodes are entangled with those of 
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FIG. 8: (Color online) The properties of the active nodes, (a) For the 
giant-cluster state shown in Fig. 4, the probability that node stays in 
the giant cluster versus the node index, (b) A segment of (a) but with 
different coupling strengths nearby the LW bifurcation point e\. (c) 
For the giant-cluster state (e = 0.952) nearby the SW bifurcation 
point, the probability that node stays in the giant cluster versus the 
node index, (d) A segment of (c) under different coupling strengths 
nearby the SW bifurcation point £2- 



We go on to investigate the affection of the coupling 
strength on the locations of the active nodes. In Fig. 8(b) we 
fix the network topology and compare the node activities un- 
der different coupling strengths nearby the bifurcation point 
e%. It is found that, despite of the changes in p u \, the lo- 
cations of the active nodes are kept unchanged. That is, the 
active nodes are always the first ones to escape from the gi- 
ant cluster whenever the network is unsynchronizable. We 
have also investigated the affection of the bifurcation type on 
the locations of the active nodes. By choosing the coupling 
strength nearby the SW bifurcation e = 0.952 > £2, we plot 
in Fig. 8(c) the node attaching probability p U 2 as a function 
of the node index i. An interesting finding is that, comparing 
to the situation of LW bifurcation [Fig. 8(a)], the locations 
of the active nodes have been totally changed in Fig. (c). In 
Tab. 2 we list the detail information about the 5 most active 
nodes in Fig. 8(c), again their locations can not be predicted 
by the node degree or betweenness. Similar to the LW bifur- 
cation, the locations of the active nodes are also independent 
to the coupling strength at the SW bifurcation, as shown in 
Fig. 8(d). 



the stable nodes. Noticing that in the BA growth model node 
of higher index in general assume the smaller degree, the ob- 
servation of Fig. 8(a) therefore indicates the independence of 
the node degree to the node stability, or the inaccuracy of us- 
ing degree to characterize the node activity. Specifically, in 
Fig. 8(a) the 5 most unstable nodes, by a descending order of 
Pui, are those of degrees k = 47, 36, 26, 10, 4, respectively. 
Except the one of k = 4, all the other nodes have higher de- 
grees. Another well-known topological property of complex 
network is the node betweenness, which counts the number 
of shortest pathes that pass through each node and actually 
evaluates the node importance from the global-network point 
of view. This global-network property, however, is also inca- 
pable to characterize the active nodes. In Tab. 1 we list the 
detail information about the 5 most active nodes in Fig. 8(a), 
where the inaccuracy of node degree or node betweenness in 
characterizing the active nodes are summarized. 



TABLE I: For the attaching probability pi plotted in Fig. 8(a), we 
list the 5 most unstable nodes and try to characterize them by a set of 
topological quantities including the node index i, the attaching prob- 
ability pi, the stability rank pi .rank, the node degree ki, the degree 
rank fc;_rank, the node betweenness Bi, and the betweenness rank 
Si .rank. 



Node index i 


Pi 


Pi .rank 


k' 1 


ki .rank 


B, 


Bi .rank 


615 


0.72797 


1 


5 


39^537 


1301 


280 


762 


0.74424 


2 


5 


39^537 


1375 


356 


680 


0.75416 


3 


4 


1^338 


1254 


680 


372 


0.7591 


4 


6 


538^645 


1440 


406 


938 


0.75972 


5 


4 


1^338 


1215 


159 



TABLE II: Similar to Tab. I but for the attaching probability pi plot- 
ted in Fig. 8(c). Comparing to Tab. I, one important observation is 
the changed locations of the active nodes due to the changed bifurca- 
tion type. 



Node index i 


Pi 


Pi .rank 




fci_rank 


B t 


Bi .rank 


43 


0.78196 


1 


9 


779^813 


2847 


813 


35 


0.78969 


2 


18 


936^940 


8513 


953 


714 


0.795 


3 


4 


1^338 


1215 


158 


130 


0.79652 


4 


13 


886^901 


4200 


886 


154 


0.19944 


5 


10 


814^846 


2898 


815 



Previous studies about network synchronization have 
shown that, while individually it is difficult to predict the dy- 
namical behavior of each node, the average performance of 
an ensemble of nodes of the same network properties do have 
some reliable characters. For instance, it has been shown that 
in complex networks the high-degree nodes are on average 
more synchronizable than the low-degree ones [19]. Regard- 
ing to the problem of node activities, it is natural to ask the 
similar question: are the high-degree nodes more synchro- 
nized than the low-degree nodes? In Fig. 9 we plot the av- 
erage attaching probability (p u i)k as a function of degree k. 
Still, we can not find a clear dependence of (pui)^ on k. 

VII. DISCUSSIONS AND CONCLUSION 

It is worthy of note that our studies of active nodes are only 
focused on the giant-cluster state, and the purpose is to under- 
stand their dynamics and reveal their properties. By ensemble 
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networks. In Ref. ^ the authors have discussed the transient 
process of global synchronization in complex networks, but 
their study are concentrating on the synchronizable state in 
which, during the course of system evolution, small clusters 
integrate into larger clusters monotonically and finally reach 
the synchronization state. After that, the system will always 
stay on the synchronization state. Our works are also different 
to the studies of Refs. [ 10, 20]. Similar to our works, in these 
studies the authors also consider the problem of pattern for- 
mation in unsynchronizable networks, but their interests are 
focused on the stationary pattern of the system. That is, the 
size and contents of the clusters do not change with time. In 
contrast, in our studies both the size and contents of the clus- 
ters are variable. 



FIG. 9: (Color online) The average attaching probability (pui) k as 
a function of node degree k. On average, the 5 most unstable nodes 
are those of degrees k — 47, 36, 26, 10, 4. Still, we can not find a 
clear dependence between (pui) k and k. 

average, we may able to improve our prediction of the active 
nodes, say for example the dependence of (p u i)k on ^ m Fig- 
9 may be smoothed if we average the results over a large num- 
ber of network realizations. Such an improvement, however, 
comes at the cost of the decreased prediction accuracy due to 
the increased candidates. Taking Fig. 9 as an example, al- 
though it is noticed that nodes of fc = 4 in general are more 
active than those of other degrees, only one of them is listed 
as the 5 most unstable nodes [Tab. 1]. In specific, among the 
total number of 338 nodes which have degree fc = 4, most 
of them are tightly attracted to the giant cluster (90 percent of 
them have attaching probabilities p u i > 0.95). Therefore, in 
terms of precise predication, the average method is infeasible 
in practice. 

Beside node degree and betweenness, we have also checked 
the dependence of the property of node activity to some other 
well-known network properties such as the clustering coeffi- 
cient, the modularity, and the assortativity. However, none of 
them is suitable to characterize the active nodes, their perfor- 
mance is very similar to that of the node degree described in 
Tab. 1 and Tab. 2. Our study thus suggests that, to give a pre- 
cise prediction to the active nodes, we may need to develop 
some new quantities. 

Despite of the amount of studies carried on network syn- 
chronization, to the best of our knowledge, we are the first to 
study the nonstationary pattern in unsynchronizable complex 



In summary, we have reported and investigated a kind of 
new phenomena in network synchronization: the nonstation- 
ary pattern. That is, the final state of the network settles nei- 
ther to the synchronization state nor to any stationary state 
of fixed pattern, the system is travelling among all the possi- 
ble patterns in an intermittent fashion (the pattern can be of 
any configuration, but its probability of showing up is pattern- 
dependent). We attribute this nonstationarity to the asymmet- 
ric topology of the complex networks, and its dynamical ori- 
gin can be understood from the property of the finite-time 
Lyapunov exponent associated to the desynchronized mode. 
Two types of synchronization formats, the complete synchro- 
nization and the temporal phase synchronization, have been 
employed to detect the nonstationary dynamics. For coupling 
strength immediately out of the stable region, the pattern evo- 
lution is characterized by the process of on-off intermittency 
and the existence of the giant-cluster; while if the coupling 
strength is far away from the bifurcation points, the pattern 
evolution is signatured by the random interactivities among 
the number of small clusters. A remarkable finding is that, 
in the giant-cluster state the locations of the active nodes are 
independent of the coupling strength but are sensitive to the 
bifurcation types. The active nodes, however, can not be char- 
acterized by the currently known network properties, further 
investigations about their identification are necessary. While 
we are hoping our studies about nonstationary pattern could 
give some new understandings to the dynamics of coupled 
complex systems, we also hope that our findings about unsyn- 
chronizable networks could be used to some practical prob- 
lems where system maintains their normal functions only un- 
der the unsynchronizable states, for example the problem of 
epileptic seizers 1I2TI1 . 
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